Mixed Model

Nick Williams
November 21, 2018

Random intercept models

no_sp_model <- organ_sp %>% 
  lmer(eligible_population_enrolled ~ total_days + opo + (1 | county), data = ., 
       REML = FALSE)

no_sp_model %>% 
  broom::tidy() %>% 
  filter(term == "total_days") %>% 
  mutate(estimate_year = estimate * 365, 
         se_year = std.error * 365, 
         ci_low = estimate_year - 1.96*se_year, 
         ci_upper = estimate_year + 1.96*se_year) %>% 
  mutate(term = ifelse(term == "total_days", "Time (year)", term)) %>% 
  select(term, estimate_year, se_year, ci_low, ci_upper) %>% 
  rename("Term" = term, 
         "$\\beta$" = estimate_year, 
         "Standard error" = se_year, 
         "Lower bound" = ci_low, 
         "Upper bound" = ci_upper) %>% 
  knitr::kable(digits = 3, escape = F) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))
Term *β* Standard error Lower bound Upper bound
Time (year) 2.675 0.008 2.66 2.69
organ_sp %>% 
  ggplot(aes(x = total_days, y = eligible_population_enrolled)) + 
  geom_jitter()

add_predictions(organ_sp, no_sp_model) %>% 
  ggplot(aes(x = total_days, y = pred, color = opo, type = county)) + 
  geom_line() + 
  geom_jitter(aes(x = total_days, y = eligible_population_enrolled), color = "black", alpha = 0.1) +
  viridis::scale_color_viridis(discrete = TRUE) + 
  theme(legend.position = "bottom")

sp_12_model <- organ_sp %>% 
  lmer(eligible_population_enrolled ~ total_days + year_sp_2012 + (1 | county), data = .)

sp_12_model %>% 
  broom::tidy() %>% 
  filter(term %in% c("total_days", "year_sp_2012")) %>% 
  mutate(estimate_year = estimate * 365, 
         se_year = std.error * 365, 
         ci_low = estimate_year - 1.96*se_year, 
         ci_upper = estimate_year + 1.96*se_year) %>% 
  mutate(term = ifelse(term == "total_days", "Time (year)", term), 
         term = ifelse(term == "year_sp_2012", "2012 spline", term)) %>% 
  select(term, estimate_year, se_year, ci_low, ci_upper) %>% 
  rename("Term" = term, 
         "$\\beta$" = estimate_year, 
         "Standard error" = se_year, 
         "Lower bound" = ci_low, 
         "Upper bound" = ci_upper) %>% 
  knitr::kable(digits = 3) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))
Term *β* Standard error Lower bound Upper bound
Time (year) 2.883 0.022 2.840 2.926
2012 spline -0.326 0.032 -0.389 -0.262
add_predictions(organ_sp, sp_12_model) %>% 
  ggplot(aes(x = total_days, y = pred, color = opo, type = county)) + 
  geom_line() +
  viridis::scale_color_viridis(discrete = TRUE) + 
  theme(legend.position = "bottom")

sp_16_model <- organ_sp %>% 
  lmer(eligible_population_enrolled ~ total_days + year_sp_2016 + opo + (1 | county), data = .) 

sp_16_model %>% 
  broom::tidy() %>% 
  filter(term %in% c("total_days", "year_sp_2016")) %>% 
  mutate(estimate_year = estimate * 365, 
         se_year = std.error * 365, 
         ci_low = estimate_year - 1.96*se_year, 
         ci_upper = estimate_year + 1.96*se_year) %>% 
  mutate(term = ifelse(term == "total_days", "Time (year)", term), 
         term = ifelse(term == "year_sp_2016", "2016 spline", term)) %>% 
  select(term, estimate_year, se_year, ci_low, ci_upper) %>% 
  rename("Term" = term, 
         "$\\beta$" = estimate_year, 
         "Standard error" = se_year, 
         "Lower bound" = ci_low, 
         "Upper bound" = ci_upper) %>% 
  knitr::kable(digits = 3) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))
Term *β* Standard error Lower bound Upper bound
Time (year) 2.563 0.01 2.542 2.583
2016 spline 0.789 0.05 0.692 0.886
add_predictions(organ_sp, sp_16_model) %>% 
  ggplot(aes(x = total_days, y = pred, color = opo, type = county)) + 
  geom_line() + 
  viridis::scale_color_viridis(discrete = TRUE) + 
  theme(legend.position = "bottom")

sp_17_model <- organ_sp %>% 
  lmer(eligible_population_enrolled ~ total_days + year_sp_2017 + opo + (1 | county), data = ., 
       REML = FALSE)

sp_17_model %>% 
  broom::tidy() %>% 
  filter(term %in% c("total_days", "year_sp_2017")) %>% 
  mutate(estimate_year = estimate * 365, 
         se_year = std.error * 365, 
         ci_low = estimate_year - 1.96*se_year, 
         ci_upper = estimate_year + 1.96*se_year) %>% 
  mutate(term = ifelse(term == "total_days", "Time (year)", term), 
         term = ifelse(term == "year_sp_2017", "2017 spline", term)) %>% 
  select(term, estimate_year, se_year, ci_low, ci_upper) %>% 
  rename("Term" = term, 
         "$\\beta$" = estimate_year, 
         "Standard error" = se_year, 
         "Lower bound" = ci_low, 
         "Upper bound" = ci_upper) %>% 
  knitr::kable(digits = 3) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))
Term *β* Standard error Lower bound Upper bound
Time (year) 2.561 0.009 2.543 2.579
2017 spline 1.700 0.078 1.548 1.852
add_predictions(organ_sp, sp_17_model) %>% 
  ggplot(aes(x = total_days, y = pred, color = opo, type = county)) + 
  geom_line() +
  viridis::scale_color_viridis(discrete = TRUE) + 
  theme(legend.position = "bottom")

all_sp_model <- organ_sp %>% 
  lmer(eligible_population_enrolled ~ total_days + year_sp_2012 + year_sp_2016 + year_sp_2017 + opo + 
         (1 | county), data = ., REML = FALSE)

all_sp_model %>% 
  broom::tidy() %>% 
  filter(term %in% c("total_days", "year_sp_2012", "year_sp_2016", "year_sp_2017")) %>% 
  mutate(estimate_year = estimate * 365, 
         se_year = std.error * 365, 
         ci_low = estimate_year - 1.96*se_year, 
         ci_upper = estimate_year + 1.96*se_year) %>% 
  mutate(term = ifelse(term == "total_days", "Time (year)", term), 
         term = ifelse(term == "year_sp_2012", "2012 spline", term), 
         term = ifelse(term == "year_sp_2016", "2016 spline", term), 
         term = ifelse(term == "year_sp_2017", "2017 spline", term)) %>% 
  select(term, estimate_year, se_year, ci_low, ci_upper) %>% 
  rename("Term" = term, 
         "$\\beta$" = estimate_year, 
         "Standard error" = se_year, 
         "Lower bound" = ci_low, 
         "Upper bound" = ci_upper) %>% 
  knitr::kable(digits = 3) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped")) %>% 
  footnote(general = "Splines correspond with dates in which 3 major policy changes were inacted. \nEstimates are adjusted for organ donor organization")
Term *β* Standard error Lower bound Upper bound
Time (year) 3.088 0.022 3.045 3.131
2012 spline -0.982 0.041 -1.063 -0.902
2016 spline 0.333 0.150 0.039 0.627
2017 spline 2.465 0.213 2.048 2.883
Note:
Splines correspond with dates in which 3 major policy changes were inacted.
Estimates are adjusted for organ donor organization
add_predictions(organ_sp, all_sp_model) %>% 
  ggplot(aes(x = total_days, y = pred, color = opo, type = county)) + 
  geom_line() +
  viridis::scale_color_viridis(discrete = TRUE) + 
  geom_vline(xintercept = 1491, linetype = "dashed") + 
  geom_vline(xintercept = 2799, linetype = "dashed") + 
  geom_vline(xintercept = 3088, linetype = "dashed") + 
  annotate(geom = "text", x = 1720, y = 5, label = "2012") +
  annotate(geom = "text", x = 2600, y = 5, label = "2016") +
  annotate(geom = "text", x = 3300, y = 9, label = "2017") +
  theme(legend.position = "bottom", 
        legend.title = element_blank()) + 
  guides(color = guide_legend(ncol = 2)) + 
  labs(title = "Proportion of organ donors in NY counties", 
       x = "Days since start of data collection",
       y = "% enrolled as donors")

Comparing nested models

mods <- list(no_sp_model, sp_12_model, sp_16_model, sp_17_model, all_sp_model)

names(mods) <- c("no_sp_model", "sp_12_model", "sp_16_model", "sp_17_model", "all_sp_model")

rmse_fun <- function(model) {
  rmse(model, organ_sp)
}

purrr::map(mods, rmse_fun) %>% 
  as_tibble() %>% 
  gather(key = "model", 
         value = "rmse") %>% 
  mutate(model = ifelse(model == "no_sp_model", "No spline", model), 
         model = ifelse(model == "sp_12_model", "2012 spline", model), 
         model = ifelse(model == "sp_16_model", "2016 spline", model), 
         model = ifelse(model == "sp_17_model", "2017 spline", model), 
         model = ifelse(model == "all_sp_model", "All splines", model)) %>% 
  rename("Model" = model, 
         "RMSE" = rmse) %>% 
  knitr::kable(digits = 3) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))
Model RMSE
No spline 1.915
2012 spline 1.902
2016 spline 1.883
2017 spline 1.855
All splines 1.770
options(knitr.kable.NA = '')

anova(no_sp_model, sp_17_model) %>% 
  broom::tidy() %>% 
  select(term, AIC, logLik, statistic, Chi.Df, p.value) %>% 
  mutate(term = ifelse(term == "no_sp_model", "No spline", term), 
         term = ifelse(term == "sp_17_model", "2017 spline", term), 
         p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>% 
  knitr::kable(col.names = c("Model", "AIC", "Log-like", "$\\chi^2$", "df", "p-value"), 
               escape = F) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))
Model AIC Log-like *χ*2 df p-value
No spline 30768.94 -15377.47
2017 spline 30306.86 -15145.43 464.084 1 < 0.001
anova(no_sp_model, all_sp_model) %>% 
  broom::tidy() %>% 
  select(term, AIC, logLik, statistic, Chi.Df, p.value) %>% 
  mutate(term = ifelse(term == "no_sp_model", "No spline", term), 
         term = ifelse(term == "all_sp_model", "All splines", term), 
         p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  knitr::kable(col.names = c("Model", "AIC", "Log-like", "$\\chi^2$", "df", "p-value"), 
               escape = F) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))
Model AIC Log-like *χ*2 df p-value
No spline 30768.94 -15377.47
All splines 29630.95 -14805.47 1143.996 3 < 0.001
anova(sp_17_model, all_sp_model) %>% 
  broom::tidy() %>% 
  select(term, AIC, logLik, statistic, Chi.Df, p.value) %>% 
  mutate(term = ifelse(term == "sp_17_model", "2017 Spline", term), 
         term = ifelse(term == "all_sp_model", "All splines", term), 
         p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  knitr::kable(col.names = c("Model", "AIC", "Log-like", "$\\chi^2$", "df", "p-value"), 
               escape = F) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))
Model AIC Log-like *χ*2 df p-value
2017 Spline 30306.86 -15145.43
All splines 29630.95 -14805.47 679.912 2 < 0.001

Random intercept and random slope models

# convergence issues with random slope, most likely not needed

# rs <- organ_sp %>% 
#   lmer(eligible_population_enrolled ~ total_days + (1 + total_days | county), data = .) 
# 
# coef(rs)$county %>% 
#   rownames_to_column(var = "county") %>% 
#   mutate(year = total_days * 365)